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A multi-grid, flux-difference-split, finite- volume code, VULCAN, is presented for 
solving the elliptic and parabolized form of the equations governing three-dimensional, 
turbulent, calorically perfect and non-equilibrium chemically reacting flows. The space 
marching algorithms developed to improve convergence rate and or reduce computa- 
tional cost are emphasized. The algorithms presented are extensions to the class of 
implicit pseudo-time iterative, upwind space-marching schemes. A full approximate stor- 
age, full multi-grid scheme is also described which is used to accelerate the convergence 
of a Gauss-Seidel relaxation method. The multi-grid algorithm is shown to significantly 
improve convergence on high aspect ratio grids. 


Introduction 

Scramjet engines and nozzles employed by hy- 
personic vehicles operate in a viscous super- 
sonic/hypersonic flow regime. Such conditions allow 
the use of space-marching solution algorithms without 
adversely affecting the fidelity of the calculation when 
the flow is steady and unseparated in the marching 
direction. Numerous algorithms have been developed 
for solving the Euler, Navier-Stokes, and Parabolized 
Navier-Stokes (PNS) equations (e.g. Refs. 1-7) for a 
thermally and calorically perfect gas. Extensions of 
space marching algorithms have also been made to 
handle both equilibrium and non-equilibrium effects 
(e.g. Refs. 8-15). In addition, the PNS equations can 
be solved much more efficiently than the Navier-Stokes 
equations, which makes the PNS equations very at- 
tractive for use in design and optimization methods. 16 

Space marching methods can be broken down al- 
gorithmically into two groups, steady-state equation 
solvers and pseudo-time iterative solvers. Steady- 
state equation solvers, as the name implies, solve 
the steady-state form of the Euler and PNS equa- 
tions. Explicit 3, 5,10 > 17 and implicit 4,11,14 steady-state 
space-marching schemes have been used with consid- 
erable success. Pseudo-time iterative space-marching 
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schemes that solve the unsteady form of the Euler and 
PNS equations 2,18,19 have also been employed using 
relaxation schemes that iterate at each spatial step in 
pseudo-time to reach a steady-state solution. 

A large number of the problems of interest can in- 
clude subsonic and/or separated flow regions, which 
requires that the space marching algorithm be imple- 
mented within the broader context of a code capa- 
ble of solving the Navier-Stokes equations. VULCAN 
(Viscous Upwind aLgorithm for Complex flow ANal- 
ysis) is a cell-centered finite- volume, structured grid, 
multi-block code which solves the equations governing 
inviscid and viscous flow of a calorically perfect gas 
or of an arbitrary mixture of thermally perfect gases 
undergoing non-equilibrium chemical reactions. VUL- 
CAN allows the flow domain to be decomposed into 
regions in which the most suitable algorithm (ellip- 
tic or marching) can be utilized. The inviscid fluxes 
are computed using the MUSCL scheme 20 with either 
the approximate Riemann solver of Roe 21 or the low 
dissipation flux splitting scheme of Edwards. 22 The 
viscous fluxes can be evaluated either with or with- 
out cross derivative contributions. VULCAN solves 
elliptic flows by marching the unsteady form of gov- 
erning equations in time. Hyperbolic flows (e.g. the 
Euler equations in supersonic flow) and parabolic flows 
(e.g. the parabolized Navier-Stokes equations) are 
solved by space marching while iterating the unsteady 
equations in pseudo-time to a steady state solution. 
The pseudo-time iterative approach was chosen be- 
cause it allows the time marching solution algorithms 
and convergence acceleration techniques developed to 
solve the elliptic Euler equation and FNS equations 
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to be used to reduce the computational cost of the 
space marching scheme. The explicit PNS algorithm 
described in a previous paper , 23 while extremely ef- 
ficient for inviscid and moderate Reynolds number 
flows, was found to be too computationally expensive 
when solving high Reynolds number flows due to the 
large number of space marching steps required by sta- 
bility constraints. The code structure of VULCAN 
was designed to allow the elliptic and space march- 
ing schemes to share numerical algorithms developed 
to solve and/or accelerate the convergence of the gov- 
erning equations to steady state. The code has in- 
corporated within it four time marching schemes: a 
multi-stage Runge-Kutta scheme, a diagonalized ap- 
proximate factorization scheme, a block approximate 
factorization scheme, and a diagonally dominant alter- 
nating direction implicit scheme with a factorization 
error correction sub-iteration. Each of these schemes 
can be used as the smoothing algorithm in a full ap- 
proximate storage (FAS) full multi-grid (FMG) scheme 
and thereby accelerate the convergence to steady state. 


Governing Equations 

The unsteady three-dimensional form of the Favre 
averaged Navier-Stokes and species transport equa- 
tions can be written in nonorthogonal curvilinear form 
as 


Qt + (E - E„) e + (F - + (G - G„) ( = S (1) 


where the subscripts £,77, and £ denote partial differen- 
tiation with respect to the nonorthogonal curvilinear 
coordinates and the subscript t denotes partial dif- 
ferentiation with respect to time. The contravariant 
components of the curvilinear nonorthogonal coordi- 
nate system are 


averaged variables (indicated as <f>) as 


pfi 


UJi 

PfN cs 


w N cs 

~pu 


0 

~pv 


0 

~pw 

s = 

0 

pE 


0 

~pk 


Sa 

pi 


Si 

pg 


s- g 

pf 


. Sf _ 


pfiU 


p/n cs U 
puu + i x p 
_ |V£| pvU + iyp 

~ J pwU + i z p 

pHU 
~pkU 

piu 

pgu 

pTU 

T pfiV 


F 


|V>?I 

J 


~P In cs V 
puV + Tj x p 

PVV + fjyf 

~pwV + TJzP 
pHV 
pkV 
plV 
pgV 
-plV 


pfiW 


( 4 ) 


( 5 ) 


( 6 ) 


£ = £(*, y,z ) ; v = v( x > y , z ) ; C = C(u y , z ) (2) 


and 


J = 


d{£,,r],C) 

d{x,y,z) 


( 3 ) 


where £ a , rj a , and G (a = x, y, or z) are the compo- 
nents of the cell face area normal vector, while £ a , £ a , 
and G are the components of the cell face unit normal. 

The Q, E, F, G, and S vectors are defined using 
Favre averaged variables (indicated as <f>) and Reynolds 
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where ~p is the mixture density; fi is the mass fraction 
and uii is the species production rate of the i th chemical 
species; N cs is the total number of chemical species; 
p is the static pressure; u, v and w are the Cartesian 
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velocity components; and E and H are the total energy 
and total enthalpy respectively. The variables k and l 
are the turbulent kinetic energy and a turbulent length 
scale equation variable (either wort). The turbulence 
chemistry interaction variables g and T are the energy 
variance e" e" , and the sum of the species variances 
f'i fi ■ The contravariant velocities are defined 
as 


u = iu + i y v + i z E 

V = r) x u + fjyl + rj z w (8) 

W = ( X U + (yV + ( Z W 


The viscous flux vectors E^ , , and G„ are written 
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The components of the species molecular diffusion ve- 
locity of the i th specie in the j th direction, Udjifi, are 
expressed using a Fickian model and their turbulent 
counterparts, u' d . f" , are closed using a gradient dif- 
fusion model as 
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where p and p t are the molecular and turbulent vis- 
cosity and Sc and Sc t are the laminar and turbulent 
Schmidt numbers. The Cartesian components of the 
viscous energy flux are 
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where (Tk is the turbulent Prandtl number for the diffu- 
sion of turbulent kinetic energy. The Cartesian viscous 
stresses are defined as 
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and the total heat fluxes are given as 
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where A is the gas mixture molecular thermal conduc- 
tivity and Pr t is the turbulent Prandtl number. The 
turbulence modeling equation diffusion terms D ^ , 

Dz and D z , where <f> is either k, c, uj , g, or T, 

9y 9 z 

are modeled as 
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Turbulence Models 

The Reynolds stress term uj") model is depen- 
dent on the type of turbulence model selected. Cur- 
rently, only two equation models are available in VUL- 
CAN. However, the two equation model implemen- 
tation allows flexibility when modeling the Reynolds 
stresses. If one of the eddy viscosity based two equa- 
tion turbulence models is selected, the Boussinesq ap- 
proximation is made where the Reynolds stresses are 
modeled as 
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If an explicit algebraic Reynolds stress two equa- 
tion turbulence model is selected, then the Reynolds 
stresses are modeled as 
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and uj = c/k. The definitions for p * , p** , oq, and oq 
can be found in Abid et al. 25 

The two equation turbulence models in VULCAN 
can be categorized as k — c based models, Wilcox’s 28 
k — uj based models and Menter’s 24 k — uj models. The 
k—c based models can be solved with either the Boussi- 
nesq model or an explicit algebraic Reynolds stress 
model. 25-27 The Wilcox k — uj based models are avail- 
able in three forms; a high Reynolds number model, 28 
a low Reynolds number model, 28 and an explicit alge- 
braic Reynolds stress model. 27 The k — uj models of 


Menter 24 are available in two forms; the baseline model 
and the Shear Stress Transport (SST) model. A cor- 
rection to each model for compressibility effects has 
also been implemented using the compressible dissipa- 
tion model of Sarkar et al. 29 for the k — c based models 
and the compressibility model suggested by Wilcox, 28 
which incorporates Sarkar’s model as well as Zeman’s 
lagging function, 30 for all of the k — uj based models 
(including Menter’s). 

Near Wall Turbulence Models 

The near wall behavior of the k — c two equation 
models is controlled through the introduction of the 
low Reynolds number modifications of Abid. 31 The 
near wall behavior of the k — uj two equation models 
is treated in either a high Reynolds number manner 
by integrating the high Reynolds number form of the 
equations to the wall without any special wall treat- 
ment or by integrating the low Reynolds form of the 
Wilcox k — uj two equation model 28 to the wall. 

If integration to the wall is not feasible for a 
given problem, the wall matching function approach of 
Wilcox 28 can be used. This approach models the wall 
shear stress and heat transfer by enforcing a compress- 
ible law of the wall which includes additional terms 
that model streamwise pressure gradients. Wilcox 32 
showed that inclusion of the pressure gradient terms 
in the wall matching function dramatically improved 
predictions of shear stress. The wall matching func- 
tion has been successfully used on grids where the y + 
(= -y/Tw/p-, t w is the wall shear stress) varied between 
0.5 and 150.0 and produced wall skin friction and heat 
transfer predictions that agreed well with results ob- 
tained on grids that were integrated to the wall for 
y + < 1.0. This is a very useful characteristic which sig- 
nificantly reduces flow solution sensitivity to the grid 
and greatly reduces the required amount of a priori 
knowledge about the range of y + when generating the 
computational grid. 

Pressure and Thermodynamic Models 

VULCAN was written to model the working gas ei- 
ther as a single component calorically perfect gas or 
as an arbitrary mixture of thermally perfect gases. 
For a single component calorically perfect gas, the 
solution vector of conserved variables in Eq. (4) is con- 
tracted to contain only a single continuity equation 
J> (with uj = 0), three momentum equations ~pu, J>v, 
and ~pw, the energy equation ~pE and the turbulence 
equations for J>k and ~pl. The pressure, static enthalpy, 
total enthalpy, and total energy for a single component 
calorically perfect gas are 

p = pRT (21) 

h = (J2L) f (22) 

H = h + ^ (u 2 + v 2 + w 2 ) + k (23) 
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p 

E = H — - 
P 

(24) 

where R is the gas constant and 7 is the ratio of specific 
heats. 

The pressure, species enthalpy, static enthalpy, total 
enthalpy and total energy for an arbitrary mixture of 
thermally perfect gases are modeled as 

i = 1 
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where 1Z U is the universal gas constant, Mi and Ri 
(= TZ u /Mi) are the molecular weight and gas constant 
of the i th species, respectively, and h° is the refer- 
ence enthalpy at a reference temperature, T 0 , of i th 
chemical species. The specific heat at constant pres- 
sure, enthalpy, and Gibbs free energy for each chemical 
species are modeled with a polynomial of static tem- 
perature. For example, the C Pt polynomial is modeled 
as either 33 


VULCAN. Wilke’s formula 35 is then applied to model 
the mixture molecular viscosity as 
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where 





(33) 


and Xi is the Hh chemical species mole fraction. 

The molecular thermal conductivity, A, of a single- 
component gas is modeled using the Prandtl number 
such that 

A = Gp ^ 

Pr 

The molecular thermal conductivity of a multi- 
component gas is modeled using either using 
Eqs. (32,33,34) or by using Sutherland’s formula for 
the chemical species molecular thermal conductivity 
A i such that 
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The coefficients A{ , Bi , C'i , Di , Ei , F( , and Gi are read 
from thermodynamic curve fit data base hies com- 
piled for each polynomial from the data in McBride 
et al. 33 > 34 


Molecular Transport Models 

The molecular viscosity, fi, of a single-component 
gas is modeled using Sutherland’s formula where 


H = Ho 


T 

T~o 


~ T 0 + S 
T+S 


(30) 


where fi 0 , T 0 , and S are constants that are specified 
for the gas of interest. The molecular viscosity, /q, of 
each component of a multicomponent gas is modeled 
using Sutherland’s formula 


where A op , T 0 ,i , and Si are constants that are specific 
to the chemical species of interest and are available as 
part of the thermodynamics data base supplied with 
VULCAN. Wassiljewa’s formula 36 is then applied to 
model the mixture molecular viscosity 
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Chemical Source Term Model 

The chemical source term, w,-, represents the mean 
of the net rate of production of the i th species in each 
chemical reaction. This source term can be modeled ei- 
ther with or without turbulence/chemistry interaction 
effects. The chemical source term, without turbu- 
lence/chemistry interaction effects, is modeled in the 
manner of , 39 such that 


T \ ~ T 0 j + Sj 

T 0 ,i j T + Si 


(31) 


where T op , and Si are constants that are specific 
to the chemical species of interest and are available as 
part of the thermodynamics data base supplied with 


No r 
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(3S) 

Eq. (38) is the production rate of the ith species as de- 
termined by the law of mass action where C m and C n 
are the participating species concentrations, v^ m and 
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Vj n their stoichiometric coefficients, and ICf } and ICb } 
the reaction rate coefficients of the forward and back- 
ward reactions. In addition, the constraint imposed 
by mass conservation exists such that = 0. 

The forward and backward reaction rate coefficients 
are calculated based on the Arrhenius law and the 
equilibrium rate coefficient, which are defined as 

K u =A 3 f^e XV (d^ (39) 

K,f. 

= jA- (40) 

The constants in the Arrhenius law, and £ a j, 

for a hydrogen-air reaction system and a complete de- 
scription of the method used to obtain the equilibrium 
rate coefficient can be found in Drummond and Hus- 
saini 37 or Carpenter. 40 

The chemical source term model which includes the 
effects of turbulence/chemistry interaction uses as- 
sumed probability density functions (PDF) to compute 
the mean species production rate. In this approach, 
the law of mass action is written as 



(41) 


The mean reaction rate coefficients, and JCbj, are 
computed using a truncated Taylor series expansion of 
the rate coefficient equations with respect to the mean 
and fluctuating static temperature, where the static 
temperature fluctuation T is obtained from the solu- 
tion of the energy variance, e" e" , transport equation. 
The higher order moments of T introduced by the 
Taylor series expansion are computed from either an 
assumed Gaussian or [3 PDF as described in Gaffney 
et al. 41,43 The mean forward and backward species 

production terms ])lm=i Cm™ and n^=i are com- 
puted from a multivariate [3 PDF where the sum of 

the species variance f'i'fi transport equation is 

solved and the variance is used in the manner described 
in Girimaji 42 and Gaffney et al. 43 

Spatial Discretization 

The governing equations, Eq. (1), are discretized 
with a cell centered, finite volume representation. The 
convection terms are treated in a second order ac- 
curate manner with the primitive variable extrapo- 
lation/interpolation scheme (MUSCL) of van Leer. 20 
Discontinuities such as shocks are captured by employ- 
ing flux limiters or by variants on the MUSCL scheme 
which employ differentiable limiters. 44,45 The numeri- 
cal flux at the cell interface may be constructed using 
a flux-difference-split scheme where the flux construc- 
tion is defined as 

Ej = i (E^ + E^) - i |A| (Q r -Q l ) (42) 


where I represents the cell interface, L and R repre- 
sent the left and right states, and A is due to Roe’s 

approximate Riemann solver. 21 A low dissipation flux 
splitting scheme based on the method of Edwards 22 
that addresses some of the entropy violation proper- 
ties of Roe’s scheme is also implemented. The low 
dissipation flux split flux construction is defined as 

E / = ai [ PL C+ E| - p R C'-E c R ]+W [D+ Pl + D~ Pr ] 

( 43 ) 

where E c is the convection contribution to the flux 
and E p is the pressure contribution to the flux. The 
parameters C± and D± (see Ref. 22) are functions 
of the left and right contravariant Mach number and 
dj. = (d£ + a R ) / 2 with a the frozen speed of sound. In 
both flux construction methods, the MUSCL scheme 
is used to reconstruct values of the primitive variables 
at the left and right states. The primitive variables 
used are the species mass fraction /), the density p, 
the Cartesian velocity components it, v, and w, the 
static pressure p, and the turbulence variables k, l, g, 
and T . 

The viscous fluxes are discretized using a finite- 
volume representation of a central difference opera- 
tor 46 which can be configured to evaluate the cell 
interface gradients using either an approximate gradi- 
ent construction where the cross-derivative terms are 
neglected or a full gradient construction in which all 
components of the derivative are included. 46 

Elliptic Solution Procedures 

The Navier-Stokes equations, Eq. (1), are hyperbolic 
in time and elliptic in space and, therefore, can be inte- 
grated in time to provide a steady state solution that is 
elliptic in space. VULCAN has several time integra- 
tion procedures from which to choose, depending on 
the type of problem being solved. A temporally sec- 
ond order accurate low storage, explicit, multi-stage 
Runge-Kutta scheme 47 can be used to integrate either 
to a steady state solution or in a time accurate man- 
ner. Alternately, due to the stiff nature of the govern- 
ing equations and the expense involved in computing 
the chemical source terms, implicit time integration 
schemes can be selected to integrate the equations 
to steady state. An Euler implicit time integration 
scheme is written in delta form as: 

[_l_ + ^E Q +«i,F Q + «i c G Q -S Q ]AQ = -R n (44) 
where E q = f§ - ff- , F q = f§ - , G q = 

H - If’ s q = lt> AQ = Qn+1 “ Qn ’ and Rn . 18 

the steady state residual at time n. The direct solution 
of Eq. (44) results in a large banded N x N system of 
equations (where N is the number of equations) which 
must be inverted at each time step. 

Solving Eq. (44) for even a medium sized three- 
dimensional problem is impractical due to storage 
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requirements. However, Eq. (44) may be factored 
to reduce the storage. Three approaches to the 
approximate factorization of Eq. (44) are currently 
implemented in VULCAN. The approach with the 
least factorization error and highest storage require- 
ments is a planar Gauss-Seidel relaxation which uses 
a Diagonally Dominant Alternating Direction Implicit 
(DDADI) 49 scheme to factor the cross stream plane. 
The forward sweep in i of a nonlinear planar Gauss- 
Seidel relaxation of Eq. (44) may be written 


[ JAt + + B,, + — Sq] AQi j/;- 

T AQj 

A^ACJqy^ — i T C^- — 

■R n /'o n_ ti o n o n 

and the backward sweep as 


(45) 
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(46) 
N matrices 


associated with replacing the S a of Eq. (44) with a 
first order difference in the a = £ ,rj, and £ directions. 


Each sweep is factored in the rj, £ plane as 


[M + (?(A, + C,)]AQE t = — R" 
[M + 9 (A f + C f )] AQ itjtk = MAQGj 
= Qi,j,k + AQ.-j.fc 

where 

M = tL + Be + 9 (Bj) + Bc) “ Sq 


(47) 


(48) 


[M + (Af + B e + Cf)] AQ*j /k = -R" 

[M + (A, + B„ + C„)] AQ**. k = MAQGj 
[M + A c + B c + C f ] AQ iJ>k = MAQ** tk (49) 
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M = 


.JAt 


- S c 
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Both factorized schemes above result in a tridiagonal 
system of block N xN matrices which must be inverted 
in either the two (rj, £) or the three (£, rj, £) directions. 
Note that Eq. (50) contains only the source term Ja- 
cobian. 

The schemes described above are storage and op- 
eration intensive. A third scheme commonly referred 
to as a diagonalized approximate factorization (DAF) 
scheme 50 significantly reduces the storage and num- 
ber of operations by diagonalizing the Jacobians and 
reducing the N x N block tridiagonals to a system of N 
scalar tridiagonals. The DAF scheme may be written 
as 

[I-JAtS Q ]AQl jtk = -JAm n 
Tf [I + JAtS e (Af )C - Af,„)] T - 1 AQ?*. k = AQN 
T„ [I + JAtS v (A„, c - A,,„)] T^AQ?*-** = AQ,t*. k 
Tf [I + JAtS c (A fiC - A fi „)] T ^ 1 AQqy k = AQ*** k 
Q&fc = + AQ.-j.fc (51) 


where T a and T“ 4 are matrices of the left and right 
eigenvectors, X a c are the inviscid eigenvalues, and X av 
are a diagonalized form of the viscous eigenvalues. 


Space Marched Equations 

The Navier Stokes equation, Eq. (1), can be made 
parabolic in space while remaining hyperbolic in time 
by dropping all diffusive terms in the streamwise (£) 
flux vector and by employing Vigneron’s 1 treatment of 
the streamwise pressure gradient terms. This results 
in an equation of the form 


and 9 is a parameter that determines the order of ac- 
curacy of the time derivative and also can be used 
to control the stability of the scheme. A value of 2 
is suggested in Ref. 51 and Ref. 52 and is used here 
as well. Note that the diagonal term M contains the 
diagonal contributions of J^Eq, S v Fq, J^Fq as well 
as Sq and is shared by both factors. This relaxation 
scheme could also be written as a conventional two fac- 
tor scheme 18 where the B^ and B^ terms are removed 
from M and combined with their A and C contribu- 
tions. However, this would increase factorization error 
without reducing cost or storage appreciably. 

A three factor approximate factorization method 
that has lower storage requirements may be written 


Q t + + (Uly 

where 


1 )E“” + (F - F„) + (G - G„) f =S 

(52) 

r o i 




J 


0 

ixp 

iyP 

tp 

0 

0 

o 

o 

o 

o 


( 53 ) 
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(58) 


The (uj v — l)Ej” term that appears in Eq. (52) is 

a correction to the streamwise flux difference, Ef, 
that parabolizes the streamwise pressure gradient. 
This correction term contains Vigneron’s streamwise 
pressure-gradient coefficient, cj v , which is constrained 
such that the hyperbolic-parabolic nature of the gov- 
erning equations is preserved with respect to the 
streamwise direction. The constraint equation that 
Vigneron obtained from an eigenvalue analysis is 


ui v = min 


cry 

l + ( 7 - 1 )Mf 


(54) 


M| = 


IP 

(£? + Cy + G ) ip It 


(55) 


where 7 is the frozen ratio of specific heats and <7 is a 
safety factor (usually = 0.95) included to account for 
nonlinearities that are not considered in the eigenvalue 
analysis. 1 

Eq. (52) is written in a correction form in a manner 
similar to Morrison and Korte. 19 Morrison and Korte 
found that the implementation of Vigneron’s stream- 
wise pressure gradient approximation in a cell centered 
finite volume context required special treatment to 
avoid the introduction of artificial pressure gradient 
terms which in turn cause errors that are strongly grid 
dependent. As in Morrison and Korte, 19 oj v is evalu- 
ated at the cell center by averaging the cell face metrics 
and cell face primitive variables obtained from either 
a first or second order fully upwind reconstruction. 


Space Marching Solution Procedures 

Following Newsome et ah, 18 Eq. (52) is solved by 
constructing the fluxes and source terms for a constant 
^ plane (three-dimensional) or line (two-dimensional 
and axisymmetric) of cells and then applying the for- 
ward sweep step of the planar Gauss-Seidel relaxation 
scheme presented in Eq. (45). Eq. (45) is modified 
slightly for the space marching procedure to become 


[ JAt ° — S Q ]AQ i ,j ifc + 


T C„AQ,-j +1 .fc+ 


A C AQ 


i ,j ,k — 1 


C C AQ 


i,j,k + l 


-Tt n (n n + 1 r > n + 1 n n 


(56) 


where is the Jacobian of the parabolized stream- 
wise flux, E^ + (u„ — 1)E^”. 

Extending the approach of Newsome et al. to in- 
clude source terms, the cross-stream plane is factored 
using a two factor scheme 


[M + (A, + B„ + C„)] AQ*j /k = -R n 
[M + (A c + B c + C f )] AQij :k = MAQ*j ik (57) 
= Q lj,k + AQ.-J.fc 


and 


M = 


I 

lAt 


+ B“ c - S 


Q 


However, an approach that has less factorization error 
is to factor Eq. (56) in the cross stream plane using 


DDADI 


[M + 0(A, + C,)]AQh fe = -r 
[M + 9 (A f + C f )] AQ.-j.fc = MAQV (59) 
Q ■ +fc = Q"j.fc + AQ.-j.fc 

and 

M = jAt +B T + 9 ( B >) + B t) “ S Q ( 6 °) 

The DDADI factorization is an expensive algorithm, 
especially when the number of equations is large. 
Therefore, a modified planar Gauss-Seidel relaxation 
which uses the diagonalized approximate factorization 
scheme is proposed that has lower storage and oper- 
ation requirements for use on moderate aspect ratio 
grids. This scheme will converge poorly when the 
streamwise grid aspect ratio becomes large; however, 
by combining the DAF factorization with wall match- 
ing functions, the wall normal grid spacing (and thus 
aspect ratio) can be relaxed and acceptable conver- 
gence behavior recovered. A planar relaxation based 
on DAF may be written as 

[i + JAt - s Q )] aqg = -jAm n 

T v [I+JAtS v (A„, c - A^)] T-'AQ?*.* = AQN fe 

T c [I + JAtS c (A CiC - A f .„)] T^ 1 AQ,-j . k = AQ**- k 

(61) 

Qqyl = Q"j.fc + AQ.-j.fc 

Either of the relaxation procedures given above is then 
used to solve each £ plane of cells beginning at the 
inflow plane and sweeping in positive £ through the 
computational domain. Each £ plane of cells is iter- 
ated in pseudo-time until convergence is reached (a 
four order reduction in the sum of the L 2 norm of the 
continuity, momentum, and energy residuals is gener- 
ally sufficient). The solution is then projected into the 
next plane of cells as an initial condition. This proce- 
dure is repeated for each plane of cells until the outflow 
plane is reached. 

Convergence Acceleration 

VULCAN has three techniques for accelerating the 
convergence of Eqs. (1) and (52). The first technique 
is local time stepping, the second technique is coarse- 
to-Hne grid sequencing, and the third technique is the 
Full Approximate Storage (FAS) multi-grid scheme of 
Brandt. 53 
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The FAS grid transfer operators for the solution, 
residual, and coarse grid corrections are those intro- 
duced by Jameson. 54 The multi-grid strategy can be 
either a V(1,0) or W(1,0) 48 type cycle. The correc- 
tions after each prolongation are smoothed using a 
scalar coefficient implicit residual smoothing to remove 
high frequency errors introduced by the prolongation 
operator. To provide a sufficiently smooth initial con- 
dition for the fine grid, the Full Multi-grid (FMG) 
method can be invoked. The use of FMG is similar 
to coarse-to-hne grid sequencing with the exception 
that multi-grid cycles are performed on each level of 
the coarse-to-hne grid sequence. 

The multi-grid scheme also has several features that 
are not commonly employed. The governing equations 
are all solved in a fully coupled fashion on all grid lev- 
els. Additionally, it was found useful for higher Mach 
number hows to construct the viscous huxes as well as 
the source terms on all coarse grids. The turbulence 
model equations are solved on all grid levels; however, 
some special treatment for wall boundary conditions 
are employed. The wall boundary condition on the 
dissipation equation is either restricted from the hne 
grid or computed on the coarse grid by evaluating the 
wall matching function on the coarse grid. 

Implementation of the FAS and FMG methods is 
identical for the hyperbolic-elliptic and the hyperbolic- 
parabolic forms of the equations with one exception. 
When solving the hyperbolic-elliptic form of the equa- 
tion, the grid is coarsened in all directions uniformly. 
However, when solving the hyperbolic-parabolic form 
of the equations, the grid is coarsened in the cross 
flow plane only. For two-dimensional and axisymmet- 
ric problems this results in grid coarsening in the rj 
direction only; while for three-dimensional problems, 
this results in coarsening in the rj and £ directions only. 

Results 

Four example problems are presented to demon- 
strate the capabilities of the multi-grid, space march- 
ing, relaxation algorithm. The first test case is in- 
viscid, calorically perfect, Mach 2 flow over a two- 
dimensional bump in a channel; the second is calor- 
ically perfect, turbulent, viscous, Mach 6 flow over 
a flat plate; the third is calorically perfect, laminar, 
Mach 3 flow over a three dimensional compression cor- 
ner; and the fourth is a two-dimensional supersonic 
turbulent diffusion flame. 

Figure 1 shows the two-dimensional computational 
grid for a Mach 2 simulation of inviscid flow over a 
bump in a channel. The grid consists of 64 cells in the 
streamwise (£) direction and 32 cells in the cross flow 
(rj) direction. The grid aspect ratio can be seen to be 
nearly one for all cells and the grid is orthogonal to the 
lower boundary. Figure 2 presents density contours 
from a space marched solution obtained using line 
Gauss-Seidel relaxation. The streamwise fluxes were 


constructed using second order, fully upwind extrap- 
olation and the cross-stream fluxes were constructed 
using MUSCL with k = 1/3 and the smooth lim- 
iter of Venkatakrishnan. 45 A grid sub-stepping feature 
available in VULCAN was employed to control oscil- 
lations near the shock. Space marching results in a 
lack of downstream information during the stream- 
wise flux construction making it impossible to employ 
flux limiters in the “classical” sense. However, VUL- 
CAN can be made to “sub-step” on the grid. Sub- 
stepping is a procedure that linearly sub-divides the 
input streamwise cells, thereby reducing the gradients 
in the streamwise direction and increasing accuracy 
while also reducing oscillations. Figure 3 shows the 
typical convergence behavior of several computational 
planes from the solution. Each plane was initialized 
with the solution from the previous plane and the CFL 
was ramped from an initial value of 100 to a final value 
of 1000 over 5 iterations. All planes converged 4 or- 
ders in the L 2 norm in 15 iterations or less. Multi-grid 
was not found to appreciably improve the already ex- 
cellent convergence rate of this problem. The problem 
was also run using the DAF based Gauss-Seidel scheme 
using a CFL ramped from 1 to 10 over 5 iterations. 
The DAF scheme required 20 iterations on average 
to reduce the L 2 four orders of magnitude. However, 
the DAF scheme had approximately the same cost per 
plane due to its lower cost per iteration. 

Figure 4 presents the Mach contours from a su- 
personic, calorically perfect (7 = 1.4) calculation of 
turbulent flow over a 0.5 meter long flat plate. The 
unit free stream Reynolds number is 2.64 x 10 7 and 
the wall temperature ratio T w /XU, = 0.65. The base- 
line k — ui turbulence model was solved to the wall on 
a computational grid of 200 cells in the £ direction and 
153 cells in the rj direction. The streamwise fluxes were 
constructed using second order, fully upwind extrap- 
olation and the cross-stream fluxes were constructed 
using MUSCL with k = 1/3 and the smooth limiter of 
Venkatakrishnan. 45 The y + of the first cell center off 
the wall was less than 0.5 for the entire plate. Figure 5 
provides the computed outflow velocity and eddy vis- 
cosity profile demonstrating that the boundary layer is 
well resolved. Figure 6 presents several typical “post- 
transition” convergence histories at a location where 
the grid aspect ratio was Pd300. The line Gauss-Seidel 
relaxation space marching scheme was used with a 
W(1,0) cycle FMG with 3 grid levels. The L 2 norm of 
the residual is reduced 12 orders of magnitude on each 
streamwise step in 32 multi-grid cycles. Although not 
shown, it is important to note that the residual of the 
k and ui equations both converged at the same rate 
and to the same level as shown here. The problem 
was run without multi-grid and was found to require 
more than 130 cycles to converge to the same level as 
shown in figure 6. 

The third case is the supersonic laminar three- 
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dimensional flow over the symmetric corner compres- 
sion flow of Korkegi and West . 55 The unit free stream 
Reynolds number for this problem is 3.07 x 10 6 and 
the wall temperature ratio is T w /XU, = 1.00. The 
computational grid, shown in Figure 7, has 64 cells 
in the £ direction and 128 cells in both the rj and £ 
directions for a total 1,048,576 cells. The grid was 
constructed such that the streamwise grid aspect ra- 
tio was 300 at the inflow plane and grew to 10,000 
at the outflow plane so that the sensitivity of the al- 
gorithm to grid aspect ratio could be examined. The 
streamwise fluxes were constructed using second order, 
fully upwind extrapolation and the cross-stream fluxes 
were constructed using MUSCL with k = 1/3 and the 
smooth limiter of Venkatakrishnan . 45 A planar Gauss- 
Seidel relaxation which used a DDADI cross-stream 
factorization was used to drive a V(1,0) cycle FMG 
with 3 grid levels. Figure 8 presents density contours 
of the corner interaction region that was computed il- 
lustrating the self similar nature of the flow. Figure 

9 presents the convergence history where the grid as- 
pect ratio was on the order of 7000. The CFL number 
was ramped from 100 to 1000 on each plane over 25 
cycles. The method did exhibit some dependence on 
grid aspect ratio. Initially, where the grid aspect ratio 
was 300, the solution required on the order of 30-40 
multi-grid cycles per plane. However, the number of 
cycles per plane gradually climbed until it stabilized 
at 50 cycles per plane. As shown in figure 9 the L 2 
norm of the residual was reduced by four orders in 
approximately 50 multi-grid cycles. When the grid as- 
pect ratio exceeded 400, the fine grid scheme without 
multi-grid required more than five times as many cy- 
cles to achieve the same level of convergence. 

The final example is a two-dimensional turbulent 
diffusion flame which was calculated by solving the 
flow over an infinitely thin splitter plate that sepa- 
rates pure hydrogen flowing at 3800 m/sec (with a 
pressure of 1 atmosphere and a temperature of 1000 
K), from a mixture of O 2 and N 2 flowing at 1200 
m/sec (with a pressure of 1 atmosphere and a tem- 
perature of 1000 K) and then solving the wake and 
shear layer that form when the fuel and oxidizer mix. 
The flow was space marched from the plate leading 
edge to the plate trailing edge solving the Menter’s 
baseline two-equation turbulence model equations to 
the adiabatic wall. The flow was then space marched 
downstream from the plate trailing edge where it de- 
velops as a wake and then as a shear layer. Figure 

10 shows the three block computational grid that was 
used. Blocks one and two consisted of 64 cells in the 
£ direction and 77 direction and block three contained 
96 cells in the £ direction and 128 cells in the 77 direc- 
tion. Blocks one and two were run non-reacting and 
block three was run reacting with a 7 species, 7 re- 
action model . 40 Figures 11 and 12 present Mach and 
water contours showing the boundary and shear layers 


as well as the flame location. The fuel and oxidizer 
can be seen to have mixed in the shear layer and the 
flame auto-ignited a short distance downstream of the 
plate trailing edge. The computation was made using 
the DAF line Gauss-Seidel scheme with a V(1,0) FMG 
cycle with 3 grid levels. Blocks one and two required 
between 70 and 100 multi-grid cycles per streamwise 
step; block three required on the order of 100 multi- 
grid cycles per streamwise step to reduce the L 2 norm 
of the residual four orders of magnitude. Typical con- 
vergence histories are shown for each of the blocks in 
figures 13, 14, and 15. 

Concluding Remarks 

The multi-grid, flux-difference split, hnite- volume 
code, VULCAN, that solves the elliptic and parab- 
olized forms of the equations governing three- 
dimensional, turbulent, calorically perfect and non- 
equilibrium chemically reacting flows was described. 
Space marching algorithms utilizing multi-grid to ac- 
celerate convergence were demonstrated. In addition, 
a Gauss-Seidel planar relaxation scheme based on di- 
agonalized approximate factorization was introduced 
and demonstrated. This algorithm was shown to re- 
duce cost and storage relative to the diagonally dom- 
inant alternating direction implicit scheme (DDADI). 
The convergence rate of the full-approximate storage 
(FAS) multi-grid scheme applied to the parabolized 
Navier-Stokes equations for calorically perfect and 
non-equilibrium chemically reacting flows was demon- 
strated and found to perform well on high aspect ratio 
grids. 
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Fig. 3 Typical convergence for space marched 
bump. 


Fig. 1 Computational grid for a bump in a chan- 
nel. 



Fig. 2 Computed density contours for Mach 2 
flow. 
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Fig. 4 Mach contours for Mach 6 turbulent flat 
plate flow. 
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Fig. 5 Outflc 
plate flow. 


Fig. 6 Conve 
flatplate. 



-3 


Mach No. 

41 2.000 




1.900 

1.800 

1.700 

1.600 

1.500 

1.400 

1.300 

1.200 

1.100 

1.000 

0.900 

0.800 

0.700 

0.600 

0.500 

0.400 

0.300 

0.200 

0.100 

0.000 


Fig. 9 Typical convergence for supersonic sym- Fig. 11 Mach contours for turbulent diffusion 
metric corner flow. flame. 
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Fig. 10 Computational grid for turbulent diffusion 
flame. 


Fig. 12 Water contours for turbulent diffusion 
flame. 
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Fig. 13 Convergence for Block 1 of turbulent dif- 
fusion flame. 


Fig. 15 Convergence for Block 3 of turbulent dif- 
fusion flame. 
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Fig. 14 Convergence for Block 2 of turbulent dif- 
fusion flame. 
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